Analysis of time-stepping methods for the monodomain model

To a large extent, the stiffness of the bidomain and monodomain models depends on the choice of the ionic model, which varies in terms of complexity and realism. In this paper, we compare and analyze a variety of time-stepping methods: explicit or semi-implicit, operator splitting, exponential, and deferred correction methods. We compare these methods for solving the bidomain model coupled with three ionic models of varying complexity and stiffness: the phenomenological Mitchell-Schaeffer model, the more realistic Beeler-Reuter model, and the stiff and very complex ten Tuscher-Noble-Noble-Panfilov (TNNP) model. For each method, we derive absolute stability criteria of the spatially discretized monodomain model and verify that the theoretical critical time-steps obtained closely match the ones in numerical experiments. We also verify that the numerical methods achieve an optimal order of convergence on the model variables and derived quantities (such as speed of the wave, depolarization time), and this in spite of the local non-differentiability of some of the ionic models. The efficiency of the different methods is also considered by comparing computational times for similar accuracy. Conclusions are drawn on the methods to be used to solve the monodomain model based on the model stiffness and complexity, measured respectively by the eigenvalues of the model’s Jacobian and the number of variables, and based on strict stability and accuracy criteria.